library(pROC)
library(PRROC)
library(caTools)
###########################################################
###GAMs: With Interaction (Figure 5, A6-A7, and Table 3)###
###########################################################
#predicted probabilities
scores<-pred.gam208$fit #data[,3]

#true binary outcomes
binary.labels <- ss.dat.f08$civconf==1 #data[,2] == 1

#threshold
threshold<-.5

#classification stats
tp <- sum( (scores > threshold) & binary.labels )
sensitivity <- tp / sum(binary.labels)
tn <- sum( (scores <= threshold) & (! binary.labels))
specificity <- tn / sum(!binary.labels)
fp<-sum( (scores > threshold) & (! binary.labels))
fn<-sum( (scores <= threshold) & binary.labels)
False.Positive.Rate<-1-specificity
False.Negative.Rate<-fn/(fn+tp)
PPV<-tp/(tp+fp)
F1<-2*((PPV*sensitivity)/(PPV+sensitivity))
F2<-(5*(PPV*sensitivity))/(4*(PPV+sensitivity))
Accuracy<-(tp+tn)/(tp+tn+fp+fn)

#print results
sensitivity
specificity
False.Positive.Rate
False.Negative.Rate
Accuracy

#Calculate AUC and Area under PR curve
colAUC(pred.gam208$fit,ss.dat.f08$civconf, plotROC=FALSE) 
pr.gam.f <- pr.curve(scores.class0 =pred.gam208$fit, weights.class0 =ss.dat.f08$civconf, curve=T)
#Plot PRC for Figure 5
jpeg("gamprc2.jpeg", width = 6, height = 6, units = 'in', res = 500)
plot(pr.gam.f)
dev.off()
#PRC AUC
pr.gam.f[[2]]
#Store all information 
Store.2.gam<-cbind(pr.curve(scores.class0 =pred.gam208$fit, weights.class0 =ss.dat.f08$civconf)[[2]],colAUC(pred.gam208$fit,ss.dat.f08$civconf, plotROC=FALSE),PPV,sensitivity,F1,F2,Accuracy)
#Table 3
Store.2.gam

###Plot ROCs (Figure A6 and A7)
#No interaction (Figure A6)
jpeg("gamroc2s.jpeg", width = 6, height = 6, units = 'in', res = 500)
roc.gam.f.208s <- plot.roc(ss.dat.f08$civconf, pred.gam208s$fit,
                           main="GAM (No Maize and Interaction)",  xlab="False positives",ylab="True positives",
                           print.auc=TRUE, print.auc.x=0.5, print.auc.y=0.1, ci=TRUE)
dev.off()

#With interaction (Figure A7)
jpeg("gamroc2.jpeg", width = 6, height = 6, units = 'in', res = 500)
roc.gam.f.208 <- plot.roc(ss.dat.f08$civconf, pred.gam208$fit,
                          main="GAM",  xlab="False positives",ylab="True positives",
                          print.auc=TRUE, print.auc.x=0.5, print.auc.y=0.1, ci=TRUE)
dev.off()

